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Abstract. Multivariate problems are typically governed by anisotropic features such as edges 
in images. A common bracket of most of the various directional representation systems which have 
been proposed to deliver sparse approximations of such features is the utilization of parabolic scaling. 
One prominent example is the shearlet system. Our objective in this paper is three-fold: We firstly 
develop a digital shearlet theory which is rationally designed in the sense that it is the digitization 
of the existing shearlet theory for continuous data. This implicates that shearlet theory provides a 
unified treatment of both the continuum and digital realm. Secondly, we analyze the utilization of 
pseudo-polar grids and the pseudo-polar Fourier transform for digital implementations of parabolic 
scaling algorithms. We derive an isometric pseudo-polar Fourier transform by careful weighting of 
the pseudo-polar grid, allowing exploitation of its adjoint for the inverse transform. This leads to a 
digital implementation of the shearlet transform; an accompanying Matlab toolbox called ShearLab 
is provided. And, thirdly, we introduce various quantitative measures for digital parabolic scaling 
algorithms in general, allowing one to tune parameters and objectively improve the implementation 
as well as compare different directional transform implementations. The usefulness of such measures 
is exemplarily demonstrated for the digital shearlet transform. 
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1. Introduction. In recent years, applied harmonic analysts have introduced 
several approaches for directional representations of image data, each one with the 
intent of efficiently representing highly anisotropic image features. Examples include 
curvelets ® E] , contourlets j5] , and shearlets [HI EI] • These proposals are inspired 
by elegant results in theoretical harmonic analysis, which study functions defined 
on the continuum plane (i.e., not digital images) and address problems of efficiently 
representing certain types of functions and operators. One set of inspiring results 
concerns the possibility of highly compressed representations of 'cartoon' images, i.e., 
functions which are piecewise smooth with singularities along smooth curves. Another 
set of results concerns the possibility of highly compressed representations of wave 
propagation operators. In 'continuum theory', anisotropic directional transforms can 
significantly outperform wavelets in important ways. 

Accordingly, one hopes that a digital implementation of such ideas would also de- 
liver performance benefits over wavelet algorithms in real-world settings. Anticipated 
applications include [5D], where missing sensors cause incomplete measurements, and 
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the problem of texture/geometry separation in image processing - for example in as- 
tronomy when images of galaxies require separated analyses of stars, filaments, and 
sheets [2U1 [TU] . 

In many cases, however, there are no publicly available implementations of such 
ideas, or the available implementations are only sketchily tested, or the available im- 
plementations are only vaguely related to the continuum transforms they are reputed 
to represent. Accordingly, we have not yet seen a serious exploration of the poten- 
tial benefit of such transforms, carefully comparing the expected benefits with those 
delivered by specific implementations. 

In this paper we aim at providing both: 

(1) A rationally designed shear let transform implementation. 

(2) A comprehensive framework for quantifying performance of directional rep- 
resentations in general. 

For (1), we developed an implementation of the fast digital shearlet transform 
(FDST) based on a digital shearlet theory which is a very natural digitization of 
the existing shearlet theory for continuous data. Other parabolic-scaling transforms, 
for example, curvelets are inherently based on operations (rotation) which translate 
awkwardly into the digital realm. In contrast, when we consider shearlets, rotations 
are replaced by shearing, which has a natural digital realization, thus enabling a 
unified treatment for the continuum and digital realm similar to wavelets. 

The framework in (2) has three benefits. First, it provides quantitative perfor- 
mance measures which can be used to tune the parameters of our implementation, 



which is publicly available at www.ShearLab.org This allows us to specify 'recom 



mended choices' for the parameters of our implementation. Second, the same 'measure 
and tune' approach may be useful to other implemcntcrs of directional transforms. 
Third, we show a way to improve the level of intellectual seriousness in applied math- 
ematics which pretends to work in image processing. We believe that widespread 
adoption of this measure and tune framework can be very valuable, since many sup- 
posedly scientific presentations are now little more than vague, numbing 'advertising' 
or 'marketing' pitches. They could instead offer quantitative comparisons between 
algorithms, and thereby be far more informative. In fact the combination of quan- 
titative evaluation with reproducible research would be particularly effective at 
producing both intellectual seriousness and rapid progress. 

1.1. Desiderata. We start by proposing the following desiderata for the fast 
digital shearlet transform FDST and its implementation: 

[Dl] Algebraic Exactness. The transform should be based on a shearlet theory for 

digital data on a pseudo-polar grid, than merely being 'somewhat close' to 

the shearlet theory for continuous data. 
[D2] Isometry of Pseudo-Polar Fourier Transform. We introduce oversampling 

and weights to obtain an isometric pseudo-polar Fourier transform, which 

allows us to use the adjoint as inverse transform. 
[D3] Tight Frame Property. The shearlet coefficients computed by the transform 

should be based on a tight frame decomposition, which ensures an isometric 

relation between the input image and the sequence of coefficients as well as 

allows us to use the adjoint as inverse transform. This property follows by 

combining [Dl] and [D2], and allows the comparison with other transforms 

in contrast to those previous two tests. 
[D4] Time-Frequency-Localization. The spatial portrait of the analyzing elements 

should 'look like' shearlets in the sense that they are sufficiently smooth as 
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well as time-localized. Localization and smoothness in frequency domain is 

ensured by definition. 
[D5] True Shear Invariance. Since the orientation-related operator of shearlets is 

in fact the shear operator, we expect to see a shearing of the input image 

mirrored in a simple shift of the transform coefficients. 
[D6] Speed. The transform should admit an algorithm of order 0(N 2 log TV) flops, 

where A^ 2 is the number of digital points of the input image. 
[D7] Geometric Exactness. The transform should preserve geometric properties 

parallel to those of the continuum theory, for example, edges should be 

mapped to edges in shearlet domain. 
[D8] Robustness. The transform should be resilient against impacts such as (hard) 

thresholded and quantized coefficients. 

1.2. Definition of the Shearlet Transform. The main idea for the construc- 
tion of the shearlet transform with discrete parameters for functions in L 2 (R 2 ) is the 
choice of a two-parameter dilation group, where one parameter ensures the multiscale 
property, whereas the second parameter provides a means to detect directions. The 
choice for a direction sensitive parameter is particularly important, since the most 
canonical choice, the rotation, would prohibit a unified treatment of the continuum 
and digital realm due to the fact that the integer grid is not invariant under rotation. 
Shearlets parameterize directions by slope rather than angles. And the shear matrix 
does preserve the structure of the integer grid, which is key to enabling an exact 
digitization of the continuum domain shearlets. 

For each a > and s G R, let A a denote the parabolic scaling matrix and S s 
denote the shear matrix of the form 

A « = {o X) and Ss = \o l)' 

respectively. To provide an equal treatment of the x- and y-axis, the frequency plane 
is split into the four cones C\\ - C22 (see Figure [LTj ) , defined by 




{(6, 6) £K 2 :£i > 1, |6/6|>1} : ^ = 21, 

{(fi,6)eR 2 :6>i, 16/61 <l} : ^ = 11, 

{(6,6)€K 2 :£i<-l, 16/61 >1} : ^ = 22, 

{(6,6) eK 2 :6<-i, 16/61 <i} : ^ = 12. 



Let now fa G L 2 (R) be a wavelet with fa G C°°(R) and supp fa C [-4, -A] U [|,4], 
and let fa G L 2 (R) be a 'bump' function satisfying fa G C°°(R) and supp fa C [—1,1]. 
We define ip € L 2 (R 2 ) by 

i£(0 = = (1.1) 

For cone C 2 i, at scale j G No :— N U {0}, orientation s = — 2 J , . . . , 2 J , and spatial 
position m € Z 2 , the associated shearlets are then defined by their Fourier transforms 

a v (0 = 2-^(SjA 4 -4) XC21 (0e 2 "^- 5sm <«> 

= 2-^^(4-^6)^2(5 + 2^) X c 21 (6e 2 ^- 3W >, (1.2) 

where rj — (j, s, m, l) index scale, orientation, position, and cone. The shearlets for 
Cn, C21, and C22 are defined likewise by symmetry, as illustrated in Figure [i~T| and 
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we denote the resulting discrete shearlet system by 



{o~ v : rj £ Nq x {-2*, ....2 j }xZ 2 x {11, 12, 21, 22}}. 



(1.3) 



The definition shows that shearlets live on anisotropic regions of width 2~ 2 i and length 
2~ J at various orientations. 

It should be mentioned that discrete shearlets - 'discrete' referring to the set of 
parameters and not to the domain - can also be defined with respect to the dilation 
matrix A 2 -j . However, in this case the odd scales have to be handled particularly care- 
fully. The attentive reader will have also observed that recently introduced compactly 
supported shearlets 22, 2E[ do not require projecting the shearlets to the respective 
cones; however, despite other advantageous properties, they do not form a tight frame 
for L 2 (R 2 ). Finally, the generating window allows in fact more freedom than (1.1 1, 
but in this paper we restrict ourselves to this (customary) choice. 



--22 



C21 




Fig. 1.1. The cones Cu - C22 o,nd the tiling of the frequency domain induced by shearlets. 

Setting C v = [j L=ll C L , we have the following theorem from |14, Thm. 3] concern- 
ing the frame properties of the discrete shearlet system. For the definition of a tight 
(sometimes called Parseval) frame, we refer to [5]. 

Theorem 1.1 ((H)). The system ([O) is a tight frame for {/ £ L 2 (R 2 ) : 
supp f C C v }. We remark that the low frequency part can be appropriately filled in 
to obtain a tight frame for L 2 (R 2 ). 

The transform associated with this system is the discrete shearlet transform, which 
for a given function / <E L 2 (M. 2 ) is defined to be the map 

N x {-2 j , . . . , 2 1 } x Z 2 x {11, 12, 21, 22} 3 r) h-> ((/, a n )) £ C. 

It is this transform, which we aim to exactly digitize. 

1.3. Ingredients of the Fast Digital Shearlet Transform (FDST). The 

shearlet transform for continuum domain data (see Figure |1.1[ ) implicitly induces a 
trapezoidal tiling of frequency space which is evidently not cartesian. By introducing 
a special set of coordinates on the continuum 2D frequency space, the discrete shearlet 
transform can be represented as a cascade of five operations: 

• Classical Fourier transformation. 

• Change of variables to pseudo-polar coordinates. 

• Weighting by a radial 'density compensation' factor. 

• Decomposition into rectangular tiles. 

• Inverse Fourier transform of each tiles. 

Surprisingly, this process admits a natural translation into the digital domain. 
The key observation is that the pseudo-polar coordinates are naturally compatible 
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with digital image processing (compare Figure [L2j ) and perfectly suited for a digiti- 
zation of the discrete shearlet transform as a comparison with the frequency tiling 
generated by continuum domain shearlets in Figure |1.1| already visually evidences. 
Fortunately, in [T] a fast pseudo-polar Fourier transform (PPFT) is already devel- 
oped. This transform evaluates the Fourier transform of an image of size N, say, on 
a pseudo-polar grid f2 of the form 57 — il 1 U 57 2 , where 

n 1 = {(-k ■ %,k) : -f < I < f , -N < k < N}, 



AT) '"J ■ 2 — *" — 2 

fl 2 = {(k, -k • § ) : -f < i < f , -N < k < N}. 



Figure 1.2 shows an illustration of the case N = 4. For an N x N image / := 




Fig. 1.2. The pseudo-polar grid H = fi 1 U fi 2 for N = 4. 

—N/2 < u, v < N/2 — 1}, the pseudo-polar Fourier transform I of / evaluated on the 
pseudo-polar grid Q is then defined to be 

N/2-1 
u,v=-N/2 

where mo > N is an integer which for the PPFT is chosen to be mo = 2N + 1 for 
computational reasons. 

The existence of PPFT suggests that we can easily and naturally get a faithful 
FDST using this algorithm. However, besides the delicateness of digitizing the con- 
tinuum domain shearlets so that they form a tight frame on the pseudo-polar grid, 
also the use of the PPFT is not at all straightforward. The PPFT as presented [T] 
is not an isometry. The main obstacle is the highly nonuniform arrangement of the 
points on the pseudo-polar grid. This intuitively suggests to downwcight points in 
regions of very high density by using weights which correspond roughly to the density 
compensation weights underlying the continuous change of variables. In fact, we will 
show that isometry is possible with sufficient radial oversampling of the pseudo-polar 
grid; however, the weights will not be derivable from simple density compensation 
arguments. 

Summarizing, the FDST of an N x N image cascades the following steps: 

1) Application of the PPFT with an oversampling factor R in radial direction. 

2) Weighting of the function values on the pseudo-polar grid by 'density-com- 
pensation-style' weights. 

3) Decomposing the pseudo-polar-indexed values by a scaled and sheared gen- 
erating window into rectangular subbands followed by application of the 2D 
iFFT to each array. 

This is an exact analogy of the discrete shearlet transform, in which the steps of 
Fourier transformation and pseudo-polar coordinate change as well as the steps of 
decomposition into rectangular tiles and the inverse Fourier transform are collapsed 
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into one step, respectively. With a careful choice of the weights and the windows, 
this transform is an isometry as we will show. Hence the inverse transform can be 
computed by merely taking the adjoint in each step. 

1.4. Performance Measurement. The above sketch does not uniquely specify 
an implementation; there is freedom in choice of weights and windows. How can we 
decide if one choice is better than another one? It seems that currently researchers 
often use overall system performance on isolated tasks, such as denoising and compres- 
sion of specific standard images like 'Lena', 'Barbara', etc. However, overall system 
performance for a system made up of a cascade of steps seems very opaque and at 
the same time very particular. It seems far better from an intellectual viewpoint to 
carefully decompose performance according to a more insightful array of tests, each 
one motivated by a particular well-understood property we are trying to obtain. 

We have developed quantitative performance measures inspired by the desiderata 
we presented in Subsection |1.1| Each performance measure produces a real value 
or a real-valued curve, thus providing a standardized framework for evaluation and, 
especially, comparison. 

1.5. Relation with Previous Work. Since the introduction of directional rep- 
resentation systems by many pioneer researchers ( [41 [51 lr?l [71 1^1 IT4] ) . various numerical 
implementations of their directional representation systems have been proposed. The 
closest ones are the curvelet, contourlet, and previous shearlet algorithms, whose main 
features we now briefly survey. 

Curvelets [3]. The discrete curvelet transform is implemented in the software 
package CurveLab, which comprises two different approaches. One is based on uncq- 
uispaced FFTs, which are used to interpolate the function in the frequency domain 
on different tiles with respect to different orientations of curvelets. The other is based 
on frequency wrapping, which wraps each subband indexed by scale and angle into 
a fixed rectangle around the origin. Both approaches can be realized efficiently in 
0(N 2 log N) flops with N being the image size. The disadvantage of this approach is 
the lack of an associated continuum domain theory. 

Contourlets [S] . The implementation of contourlets is based on a directional filter 
bank, which produces a directional frequency partitioning similar to the one generated 
by curvelets. The main advantage of this approach is that it allows a tree-structured 
filter bank implementation, in which aliasing due to subsampling is allowed to exist. 
Consequently, one can achieve great efficiency in terms of redundancy and good spatial 
localization. A drawback of this approach is that various artifacts are introduced and 
that an associated continuum domain theory is missing. 

Shearlets [T2J [55]. In [12], Easley et. al. implemented the shearlet transform 
by applying the Laplacian pyramid scheme and directional filtering successionally. 
One drawback is the deviation from the continuum domain theory. Another draw- 
back is that the associated code was not made publicly available. In contrast to this 
implementation which is based on bandlimited subband tiling - similar to the imple- 
mentation of curvelets in [3] - in |28j , Lim provided an implementation of the shearlet 
transform based on compactly supported shearlet systems (see also [22]). These com- 
pactly supported shearlets are separable and provide excellent spatial localization. 
The drawback is that they do not form a tight frame, hence, the synthesis process 
needs to be performed by iterative methods. We further wish to mention two novel 
approaches [27J and [H] for which however no implementation is yet available nor was 
their focus on deriving an exact digitization of the continuum domain transform. 
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Summarizing, all the above implementations of directional representation systems 
have their own advantages and disadvantages, one of the most common shortcomings 
being the lack of providing the unified treatment of the continuum and digital world. 
Our effort will now be put to provide a natural digitization of the shearlet theory (ban- 
dlimited shearlets) fulfilling the unified treatment requirement as well as a software 
package ShearLab quantifying performances of directional representation systems. 

1.6. Contribution of this Paper. The contributions of this paper are two- 
fold. Firstly, we introduce a fast digital shearlet transform (FDST) which is rationally 
designed based on a natural digitization of shearlet theory. Secondly, we provide a 
variety of quantitative performance measures for directional representations, which 
allow tuning and comparison of implementations. Our digital shearlet implementation 
was tuned utilizing this framework, so we can provide the user community with an 
optimized representation. 



All presented algorithms and tests are provided atjwww . ShearLab . org in the spirit 
of reproducible research 

1.7. Contents. Section [2] introduces the fast digital shearlet transform FDST 
and proves isometry. In Section [3j we then discuss two variants of an inverse digital 
shearlet transform, namely, a direct and an iterative approach. In Section [4j we 
prove several mathematical properties of the FDST such as decay properties of digital 
shearlet coefficients. The following section, Section [5j is concerned with details of 
the associated ShearLab implementation at www. ShearLab.org. The FDST is then 
analyzed in Section [7] according to the quantitative measures introduced in Section [6j 

2. FDST for Finite Data. We start by discussing the three steps in the FDST 
as described in Subsection |1.3[ which we for the convenience of the reader briefly 
repeat: 

1) Application of the PPFT with an oversampling factor R in radial direction. 

2) Weighting of the function values on the pseudo-polar grid by 'density-com- 
pensation-style' weights. 

3) Decomposing the pseudo-polar-indexed values by a scaled and sheared gen- 
erating window into rectangular subbands followed by application of the 2D 
iFFT to each array. 

We will also show that careful selection of the oversampling factor, of the weights, 
and of the windows yields an isometric transform, which enables us to compute the 
inverse shearlet transform by its adjoint (see Section [3]). 

2.1. Weighted Pseudo-Polar Fourier Transform. Given an N x N image 
/, it is well known that the Fourier transform I of / evaluated on a rectangular N x N 
grid is an isometry: 

JV/2-1 N/2-1 

£ \I(u,v)\ 2 = — Yl \H"*,"y)\ 2 - (2-1) 

u,v=-N/2 w x ,w y = -JV/2 

This is the Plancherel formula for a function defined on a finite group |21j . 

We now intend to obtain a similar formula for the Fourier transform of / evaluated 
on the pseudo-polar grid. For this, we first extend the definition of the pseudo-polar 
grid slightly by introducing an oversampling parameter R > in radial direction. 
This new grid, which we will denote in the sequel by £Ir, is defined by 



n b = (ii u n 2 
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where 



U_2k 'M '2k 
U R ' N ' R 

{ ^ If 77 



< /> < A 

2 — — 2 



A' 



flJV < fc< 



A r / 



JY 
' 2 



A' 



2 
.RjV 
2 



< k < 



RN 



}• 



(2.2) 
(2.3) 



This grid is illustrated in Figure 2.1 Notice that the 'original' pseudo-polar grid (see 




Fig. 2.1. The pseudo-polar grid tt R = Q R U Q 2 R for N = 4 and R = 4. 



Figure 1.2 1 as introduced in Subsection 1.3 is a special case of this definition when 
choosing R — 2. Also observe that the center 

C = {(0,0)} 

appears N + 1 times in £7^ as well as 0|j, and the points on the seam lines 

4 = {(-!,f) : -f <Kf-M0}, 
4 = {(!,-f):-f<Kf,M0}, 

appear in both Sl R and f2|j. Later, we will also utilize a further partitioning of the 
sets ft R and £l 2 R as 



n]} U C U ft I 2 and 



ft| = ftl 1 ucu ft 2 !, 



where 



^ = {(-f 


2£ 
AT' 


2fe \ 
ft I 


N 
2 


<(■< 


N 
2 ' 


l<fc< 






2t 
N 1 


2k\ 
R > 


N 
2 


<i< 


N 
2 ' 




1}, 




2fe 
ft 


2i\ 
N> 


N 
2 


<i< 


JV 
2 ' 


1 < k < ^f} 




^ 2 = {(f,- 


2fc 
ft 


2i\ 
N > 


N 
2 


<i< 


N 
2 ' 


-Zf <fc<- 


!}■ 



Now our goal is to choose weights w : £Ir —> K + so that, for any N x N image /, 

w(oj x ,U) y ) ■ \i(uj x ,ujy)\ 2 , (2.4) 



N/2-1 

E 

u,v=-N/2 



\iM?= E 



where here we modify the definition of the Fourier transform according to |1J and 
define it by 



l(uJ X ,UJy 



N/2-1 

E 

u,v=-N/2 



I(u,v)e m ° 



(uojx+vujy) 



(UJ X ,LOy) € Q R , 



(2.5) 
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where uiq > N. Also notice that the factor 1/7V 2 appearing in (2.1) will now be 
hidden in the weights w (uj x , uj y ) . 

We start by computing the right hand side of 12.4b: 



} J W{U) X ,U V ) ■ \I(L) x ,UJy)\ 2 



= ^ w(UJ X ,UJy) 

(LJ x ,UJy)£Q R 

= ^ W(UJ X ,U)y) 



N/2-1 

I(u,v)e m o y x v ' 

u,v = -N/2 

N/2-1 N/2-1 

u,v=-N/2 u' ,v' = -N/2 
N/2-1 

w(uJ x ,U} y )- ^ \I(U,V)\ 2 
(uj x ,ujy)£Q,pi u,v— — N/2 

N/2-1 

+ J2 iMWy) 



^{{u-u')u x -\-{v-v')<jJ y ) 



u,v,u' ,v' = -N/2 



w(uj x ,uj y )-e m o 



- ((w— u')lj x -\-(v — v')uj y ) 



Choosing / = c UljVl S(u — U\,v — vi) + c U2jV2 S(u — u 2 ,v — v 2 ) for all —N/2 < ux, Vi, 
u 2 , v 2 < N/2 — 1 and for all c„ ljt)1 , c U2tV2 € C, we can conclude that ( |2.4[ ) holds if and 
only if 

w(u: Xl u y ) . e -^ u "' +vu * ) =6(u,v), -N + l<u,v<N-l. (2.6) 



This is equivalent to the two conditions 

Y w(uj x ,oj y ) ■ cos(^(uu) x +vu)y)) = S(u,v), 

Y w(w x , u) y ) ■ sm(^(uw x + vu)y)) = 0, 



(2.7) 
(2.8) 



for all — N + 1 < u, v < N — 1. In view of the symmetry of the pseudo-polar grid, it 
is natural to impose the following symmetry conditions on the weights: 

[SI] W(UJ X ,UJ V ) = W{Uly,U X ), (UJ X ,UJy) e tt R , 

[S2] w(UJ X ,U>y) = w(~UJy,UJ X ), (U) X ,L}y) G d R , 
[S3] W(W X ,Uy) = w(-U) X ,U)y), (Uj x ,0jy) £ tl R , 

[S4] w(uj x ,ui y ) =w(u x ,-u) y ), {u) x ,u y ) € Vt R . 
In this case, (2.8| automatically holds. By the sum formula for trigonometric func- 
tions, ( |2.7[ ) is then equivalent to 

H ' [ cos (^«^)cos(^Wj / ) - sm(^wj x )sin(^VL) y )] = S(u,v) 

for all — N + 1 < u, v < N — 1. Again, by the symmetry of the weights, this is 
equivalent to 



H w (^^v) ■ l C0S (^ UUJ x) cos{^vuj y )} =S(u,v) 



(2.9) 
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for all — N + 1 < u, v < N — 1. This is a linear system of equations with RN 2 /4 + 
RN/2 + 1 unknows and (27V — l) 2 equations, wherefore, in general, we need the 
oversampling factor R to be at least 16 to enforce solvability. 

For symmetry reasons, we can now restrict our attention to one quarter of a 



cone, say Qfj. Using (2.2), (2.3 1 , and (2.9|, we then obtain the following equivalent 
condition to (2.6 1: 



RN/2 



2k y 
mnR 1 



6(u,v)=w(0,0) + i- ]T E ™(l>-§ -w)- cos ( 2 ™ 

1=0, N/2 fc=l 
N/2-1 RN/2 



cos(2tti> ■^ k s-%) 



2k _ 21 
rn a R ' N 



) (2-10) 



i=i k=i 



for all — N +1 < u,v < N — 1. Concluding, we have the following result. 

Theorem 2.1. Let N be even, let flu — f^Ufi^ be the pseudo-polar grid defined 
in (2.2) and ( |2.3[ ) ; and let W : Qr — » R + be a weight function satisfying the symmetry 
conditions [SI] - [S4]. Then 



N/2-1 

E i j mi 2 

u,v=-N/2 



w(ui x ,u>y) ■ \l(u x ,uj v )f 



ZioZds i/ and on/j/ i/ £/ie weights w(u) x ,u> y ), (ui x ,u> y ) € f2_R satisfy condition (2.10). 
Moreover, in general, R needs to be at least 16 /or smc/i weights to exist. 

2.2. Weight Functions. To avoid high complexity in the computation of the 



weights satisfying Theorem 2.1 we relax the requirement for exact isometric weight- 
ing. Instead of representing the weights as the solution of a large system of equations, 
they will be represented in terms of an undercomplete basis for functions on the 
pseudo-polar grid. More precisely, we first design basis functions Wx , . . . , w n : fin — > 
R + such that X)j=i Wji^x^y) 7^ for all (u) X) U) y ) 6 We then define the weight 



function w : Qf 



1+ to be 



En 
3=1 ' 



with c\ , . . . , c n being nonnegative 



constants. These coefficients are determined by solving (2.10) with respect to this 



weight function w using the least square method. We compute the coefficients in this 
expansion once for a given problem size; then hardwire them in the algorithm. 

2.2.1. Recommended Choices of Weights. In what follows, we present sev- 
eral designs of basis functions, each one providing nearly isometric weighting. Notice 
that, slightly abusing notation, we will use (u) x ,u} y ) and (k,£) interchangeably. The 
weighting for the three choice we recommend is displayed in Figure |2.2| 

Choice 1: Our first choice are the following seven functions Wx, . . ■ ,w 7 : 
Center, wx = l(o,o) and w 2 = l{( Ux ,u y y\k\=i}, 

Boundary: w 3 = l{{u x ,u y y\k\=NR/2,u x =w y } and w 4 = \{u} x ,u y y.\k\=NR/2,u x ^u y }, 
Seam lines: w 5 = \k\ ■ l{(u, x ,u y ):i<\k\<NR/2,u> x =u y }> w 6 = l{(w x ,w y ):|fc|=Affl/2-3,w x =w y }, 
Interior: w 7 = \k\ ■ l{(u X! u v y.i<\k\<NR/2,u xl iuy}- 



Choice 2: This is a simplified version of 'Choice 1' using the 5 functions: 
Center: Wx — 1(o,o)j 

Boundary: w 2 = l{(u xtUy y. 1*1=^/2,^=^} and w 3 = l{{ Ux ,u y ):\k\=NR/2,u x ^ v }^ 
Seam lines: w 4 = \k\ ■ l{(u x ,u v ):i<\k\<NR/2,u x =u v }, 
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Interior. w 5 = \k\ ■ l{( U:c ,oj v y.i<\k\<NR/2,w x ^u, v }- 

Choice 3: Finally, we suggest the following N/2 + 2 functions on the pseudo- 
polar grid: 
Center. w\ = l(o,o)j 

Radial Lines: w i+2 = l {{ ^^ y y 1<]k]<NR/2 Uy= t Uic}l I = 0, 1, . . . , N/2. 



\il i_J H 



Fig. 2.2. Recommended choices of weighting of the pseudo-polar grid for N = 128 and R - 
The Choice 1-3 are displayed from left to right. 



2.2.2. Comparison of Weight Functions. The patterns of the weights are 



seemingly similar in view of Figure 2.2 However, their performances can be quite dif- 
ferent depending on the chosen performance measure. We mention that the measures 
chosen below are also part of our framework of performance measures for parabolic 
scaling algorithms discussed in Section [6j 

Letting R — 8, we generate a sequence of 5 random images 7i, . . ., I5 of size 
N x N with standard normally distributed entries. We use the following measure to 
compare the performance of different weights: 



Mi, 



1 5 



\P*wPI i 



1=1 



J -t J u 



where P : I — >• / denotes the PPFT from (2.5) and w 
- denotes the 'weighting operat or' J, 
function w : ftp —> M + . Table 2.2.2 displays the performance of the weights from 
Choice 1-3 with respect to this measure. 



by abusing notation 
w ■ J, for a image J : Or — > C and a weight 



Table 2.1 
Comparison using random images 



N 


32 


64 


128 


256 


512 


Choice 1 


4.3E-3 


2.6E-3 


2.2E-3 


1.4E-3 


9.3E-4 


Choice 2 


4.2E-3 


4.0E-3 


1.8E-3 


1.5E-3 


8.8E-4 


Choice 3 


9.8E-3 


6.2E-3 


3.4E-3 


2.1E-3 


N/A 



Next, we choose the real image 'Barbara', which we denote by /, and the measure 
^ P f^THi^ *° com P are the performance of the different choices of weights. 

From Tables |2.2.2| and |2.2.2| it can be seen that the operator P*wP seems to 
converge to an identity operator as N — > 00. The data also indicates that it is 
justifiable to define weights which are linearly increasing along the radial direction. 
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Table 2.2 
Comparison using 'Barbara' 



N 32 64 128 256 512 

Choice 1 2.4E-3 1.6E-3 1.1E-3 5.2E-4 2.2E-4 

Choice 2 2.8E-3 1.2E-3 8.3E-4 3.9E-4 1.5E-4 

Choice 3 5.6E-3 2.8E-3 2.2E-3 9.1E-4 N/A 



When using iterative methods, for instance, the conjugate gradient method, for 
computing the inverse, a weight w plays the role of a preconditioner. Therefore its 
performance as such can be effectively measured by the condition number of the 
operator P*wP, i.e., cond{P*wP) — \ max (P*wP)/\ m i n (P*wP). This measure is 



displayed in Table 2.2.2|for our selected three choices of weights. Notice that, for each 



choice of a weight function with the exception of Choice 3, the condition numbers of 
P*wP are always smaller than 2. 

Table 2.3 
Comparison of cond(P*wP) 



N 


32 


64 


128 


256 


512 


Choice 1 


1.328 


1.483 


1.621 


1.726 


1.834 


Choice 2 


1.379 


1.503 


1.621 


1.731 


1.833 


Choice 3 


1.760 


1.887 


2.001 


2.104 


N/A 



2.3. Windowing. According to our discussion of the main steps of the FDST in 
Section [2j after performing a weighted pseudo-polar Fourier transform, the data has 
to be windowed using scaled and sheared versions of a generating window function. 

In this section, we now define such a set of window functions, which we will coin 
digital shearlets, and prove that these - similar to the continuum domain - form a tight 
frame for functions J : Qr — > C. We remark that this construction is a digitization of 



the continuum domain discrete shearlets introduced in P3] (compare also (1.2) and 



(1.3)). However, it is far from obvious that again a tight frame is derived, since we 
here consider a finite domain. 

For the convenience of the reader, we first briefly recall the notion of a tight frame 
in this particular situation. Let /, g : Qr — ¥ C be two functions defined on CIr. Then, 
the inner product {f,g)n R is defined to be (f,g)n R '■= J2 x en R ( x )9( x )- ^ sequence 
{ipx : D,r — > C : A G A} with A being an indexing set is a tight frame for functions 
J : fl R -> C, if 



AeA 

It then follows from basic frame theory (see [5] ) that this allows recovery of a function 
J : £Ir —> C from its coefficients ((J, <p\)n R )\eA by computing 

AeA 

Despite the danger of repeating ourselves, let us mention that our fundamental 
goal is to introduce digital shearlets as the exact digitization of continuum domain 
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shearlets. We now describe the construction step by step, which will give evidence to 
the fact that we achieved this goal. There will be one step though - when defining 
the modulation - where we have to slightly deviate from an exact digitization, and 
we will explain the reasons for this. 

We start by defining the scaling function and the generating digital shearlet. For 
this, let jh '■= — R°g4(-^/2)l , which will soon be shown to be the lowest possible scale. 
Let Wo be the Fourier transform of the Meyer scaling function such that 

supp W C [-1, 1] and W (±l) = 0, (2.11) 

and let Vo be a 'bump' function satisfying 

supp Vo C [-§, §] with V (0 = 1 for |f | < 1,£ G R. 

Then we define the scaling function <fi for the digital shearlet system to be 

0(6,6) = Wb(4-^6)Vo(4-^C a ), (6, 6) G K 2 - 

We will later restrict this function to the pseudo-polar grid. 

Let next W be the Fourier transform of the Meyer wavelet function with 

supp^C [-4,i]U[i,4] and W{±\) = W(±A) = 0, (2.12) 

as well as 

rio g4 N] 

|Wo(4- J H)| 2 + \W(4- j £)\ 2 = l for all \(\ < N, £ G R. (2.13) 

3=3 L 

We further choose V to be a 'bump' function satisfying 

supp V C [-1, 1] and V(±l) = 0, (2.14) 

and also 

|V(f - i)| 2 + |y(OI 2 + + 1)| 2 = i for all |e| <Uei. 

Notice that this implies 

\V{2 j £,-s)\ 2 = 1 for all |f | < 1, £ G R and j > 0, (2.15) 

which will become important for the analysis of frame properties. For the choice of 
Vo, Wo, V, and W in our implementation, we refer to Section [5j Then the generating 
shearlet ip is defined as 

= W(fc)V(£), (6,6) GR 2 . (2.16) 

We will now define digital shearlets on ilj^ and extend the definition to the other 
cones by symmetry. At this time, we assume R and N are both positive, even integers 
and N = 2 n ° for some integer n eN. 

For this, we first analyze the exact digitization of the coefficients of the discrete 
shearlet system from Subsection 1.2 for a function J : FIr — > C by using the shearlets 
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ip defined in (2.16). This will lead to the appropriate range of scales and to the 
support of a scaled and sheared version of the shearlet if). 

When restricting to the cone fi^ 1 , the exact digitization of the coefficients of the 
discrete shearlet system is 



w: = (w x ,LJ !) )en| ; 1 



(2.17) 



where j, s, and m are to be determined. The choice of tp leads to the coefficients 

J(Lo x ,Lj y )2-iiw{4-3Lj x )V(s + 23^)e- 2 ^ A *-> S! ' m ^ 

RN/2 N/2 

= ^K,^)2^ l M / (4^'f )V(s - 2J+ 1 ^)e- 2m ( m ' S ^ A *-^) . 

k=l i=-N/2 



The support conditions (2.121 and (2.14) of W and V, respectively, imply 



fc = 4^- 1 f+ni, m = 0,...,4 j - 1 • if 3 , 



as well as 



I = 2~ : >- 1 N(s - 1) + 



n 2 , 



n 2 = 0, 



,2--»JV, 



if we assume k and I to be positive integers. 

We next analyze the support properties in radial direction. If j < — |~log(i?/2)] , 
then k < 1, which corresponds to only one point - the origin -, and this is dealt 
with by the scaling function. Hence the lowest possible scale is jl = — [log(i2/2)] . If 
j > riog 4 AT], we have k> ^f. Hence the value W(l/4) = (cf. ( |2.12| )) is placed 
on the boundary, and thus these scales can be omitted. This implies that the highest 
possible scale is jn ■— R°g4 ^1 • Hence, the scaling parameter will be chosen to be 

j G {jL, ■ • • Jh}- 

The radial support of the windows associated with scales < j < ]h is 



fc = 4 J '- 1 f + m, ni = 0, ...A 3 



-1 15R 



(2.18) 



and the radial support of the windows associated with the scales j'l = — |~log 4 (i?/2)] 
and j H = r io S4 N 1 is 



k = 



ni, 



m = l, 



for j = jL, 



fc = 4 J '" -1 f + m, 7ii=o, 



,^-4^- 1 f, for j=j H . 



(2.19) 



We further analyze the precise support properties in angular direction. First, we 
examine the case j > 0. If s > 2 J , we have I > N/2. Hence the value V(— 1) = 
(cf. ( |2.14 1) is placed on the seam line, and these parameters can be omitted. By 
symmetry, we also obtain s > —2 J . Thus the shearing parameter will be chosen to be 

The angular support of the windows at scale j associated with shears — 2 ] < s < 2 J is 
£ = 2-1-^(3 - 1) + 7i 2 , n 2 = 0, . . . , 2-3 N, (2.20) 



SHEARLAB: A RATIONAL DESIGN OF A PARABOLIC SCALING ALGORITHM 15 



the angular support at scale j associated with the shear parameters sl = 
s H = 2 J is 



£ = 2-i- 1 N(s L -l)+n 2 , 
£ = 2-i- 1 N{s„-l)+n 2 , n 2 = 



> 2 ' 



,2~°N, £ots=sl, 



for s = sh- 



-T and 



(2.21) 



For the case j < 0, we simply let s = and £ = —N/2 + n 2 with n 2 — 0, . . . , N. 
Also, in this case, the window function W(4~ 3 lj x )V(s + 2^^-) is slightly modifed to 
be W{4~ ui x )Vq{s + 2 3 ^f-) so that the tight frame property still holds. 

These computations allow us to determine the support size of W{4~ 3 Wx)V{s + 
23 —it) in terms of pairs (k,£). In fact, the number of sampling points in radial and 
angular direction affected by a window at scale j and shear s are 



C 1 = 



zu+i R 
* i 

-1 15R , 
2 

BM. _ 47-1 

1 ~ n 



4' 



ix < i < iff) 

j = .in, 



(2.22) 



and 



CI 



-3N + 

- j a. , 

2 ' 



-2* < s < 2-? with j > 0, 
s 6 {-2^2-?} with j > 0, 

J < o, 



(2.23) 



respectively. 

Notice that the support of the continuum function £ H> I4 7 (4~ : '£;i)^(s + 2 J is of 
approximate size A J x 2 J obeying parabolic scaling. The situation is however different 
in the digital realm. Since the sampling density in angular direction does change with 
growing radius - the sampling grid becomes in fact coarser -, parabolic scaling is not 
such directly mirrored in the relation between £j and £^ s . 

Let us next carefully examine the exponential term, which can be written as 



We now adjust the exponential term as illustrated in Figure [273| which will be the only 
slight adaption we allow us to make when digitizing. The reason is to enable a direct 
application of the inverse fast Fourier transform. The necessity for this modification 
occurs because of two reasons: 

1. We cannot make the change of variables r := SjA^-jU in formula (2.17), 



which is the first step in the 'continuous' proof for tightness, due to the fact 
that the pseudo-polar grid is not invariant under the action of SjA^-j . 
2. The Fourier transform of a function defined on the pseudo-polar grid does 
not satisfy any Plancherel equation. 
Defining 8 : M. \ {0} —> M by 9(x, y) — (x, -), we let the new exponential term be 

-2 7 ri(m,(eo(S i T )- 1 )(4-^,4-^f -2-^)) = -2™(m, (4"' f ,-2^ + 1 £ )) 



This exponential term can be rewritten as 



o -2i r i{m,(4,-*%,-2i+ 1 { r )) = 2 7ri(=L+(l- S )m 2 ) -2«(m,(4"^,-2^ 1 ^)) 
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T\-l 



Fig. 2.3. Adjustment of the exponential term through the map 8 o (Sj) 



with rti and n 2 ranging over an appropriate set denned by (2.18), 

recall the < 

■CU-H2/R) 



(2.21). The reformulation - recall the definitions of and £j s in (2.22 



2l9|), and JT20) 
and |233F- 



exp 



-£ 2 2 j + 1 (l/N) 

-fix, — ^2 n 2 



nx, n 2 , 



shows that we might regard the exponential terms as characters of a suitable locally 
compact abelian group (see [21;) with annihilator identified with the rectangle 

Hi,, = {((4)-^f • rx, -(4,)-^ ■ r 2 ) : n = 0, . . . , C) - 1, r 2 = 0, . . . , t\ s - 1}. 

For the low frequency square we further require the set 

7J = {(ri,r 3 ):r 1 = -l,...,l ) r 2 = -f , . . . , f }. 

We are now ready to define digital shearlets defined on the pseudo-polar grid ftjj. 

Definition 2.2. Retaining the definitions and notations derived in this sub- 
section and Subsection 2.1, on the cone flj^ we define digital shearlets at scale j G 
{jl, ■ ■ ■ ijn}, shear s G {— 2 J , • • • , 2 J }, and spatial position m G TZj, s by 

where V J = V for j > and V 3 = Vq for j < 0, and 

1 



1 

V2 



(w^LOy) G (5* US£)\C, 

(UJ X ,LOy) e C. 



„ m on the remaining cones are defined accordingly by 



The shearlets <jf s ml a) 2 s m ,af s 
symmetry with equal indexing sets for scale j, shear s, and spatial location m. For 
lq = 1, 2 and n G 1Z, we further define the functions 



Summarizing, we call the system 
VSU = {cp% : L = 1, 2, n G ft} U {(rj )S>rn {./, a G {-2*, ■ • • , 2*'}, 



m€K jtS ,i= 11,12,21,22} 



i/ie digital shearlet system. 
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The just denned digital shearlet system forms indeed a tight frame for functions 
J : Qr — >• C, as the following results shows. 

Theorem 2.3. The digital shearlet system DSH defined in Definition \ 2. 2\ forms 
a tight frame for functions J : Or — ^ C. 

Proof. Letting J : Or C, we claim that 



(J,J) Sifl =^|(J I tf) fi J 2 + £ l(J,<.,m)n«| 



(2.24) 



which proves the result. 

We start by analyzing the first term on the RHS of (2.24). Let Lq € {1,2} and 
J c : Or -> C be defined by J c (u x ,ujy) := C(w x ,w y ) ■ J(w x ,ujy) for (ojj,^) e Or. 
Using the support conditions of </>, 



\( J 'fnh R \ 2 = | 51 J(u x ,u y )v%{u x ,w v ) 

]f\ E I E ■fcfax.Wv) • ■ e -^(«,(|^)) 



ft 



1 iV/2 



^E|E E JcK.^.^.^.e- 2 ^!.^)) 
fc=-i e=-N/2 



I ft 



(2.25) 



The choice of 7?. now allows us to use the Planchcrcl formula. Exploiting again support 
properties, we conclude from ( |2.25[ ) that 

EK J >P«> n *l a = E |CK,Ww)-^,Wy)| a -|0K,Ww)| a - 



Combining i = 1, 2 and using (2.11 1, we proved 

K^n>n*| a = E l J K-^)l 2 • l^oK)| 2 



(2.26) 



Next we study the second term on the RHS in (2.24). By symmetry, it suffices 
to consider the case t = 21. By the support conditions on W and V (see (2.12) and 
ill), 



E 



,., n /fifli 2 = E E i E j ( w *^K 2 ,w^- 

j> melius (w x ,w H )eO^ L 



E^"7 E i E Jc(u x ,u; y )-W(4-iu, x ) 



.^JTT^T^y . e -2 ff i(m,(4-^ a; ,2^)) I 2 

! 4 3 + 1 (fi/2) 2- 3 " 1 iV(s+l) 

E^ E i E E Mu*,"v) 

j,s 1 J,a ' m&n jt , k=V- 1 (R/2) i=2-i- 1 N(s-l) 

■W(4riu x ) ■ V^s + 2J^) ■ e -**i(™,(.*- i %,-*' +1 -k)) 2 . 



(2.27) 
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Similarly as before, the choice of !Zj lS does allow us to use the Plancherel formula. 
Hence (2.271 equals 



E K4<#,Jn*l a =E E Jc(^ y )-W(i-Ju } ^(s + 2^) 



Next we use (2.151 to obtain 



E E Jci^y) ■ W{A~iu x ) ■ Vi{s + 2^) 
= E \Jc(u x ,u y )\ 2 f2\W(4- j cj x )\ 2 - E \V j (s + 2^)f 

= E l^cK,^)| 2 E l^(4^)| 2 - 

(w x ,w B )eo| i :L j=jx, 



Hence the second term on the RHS in (2.24 1 equals 



JH 



E E K4"i,.,m>nJ 2 = E \J(^^v)\ 2 ■ E W 



V.)l 2 . 



(2.28) 



Finally, our claim ( |2.24[ ) follows from combining ( |2.26| , ( |2.28[ ), and ( |2.13[ ). □ 

3. Inverse FDST for Finite Data. In this section, we will analyze the inverse 
of the FDST. For this, let P and w denote the operators for the pseudo-polar Fourier 
transform and the weighting defined as before. By slight abuse of notation, we will 
further let W denote the windowing with respect to the digital shearlets, i.e., for 
J : fifl -> C, 

WJ :={{J,<p%)n R : i ,n}U{(J,a L js m ) nR :t,j,s,m} 



with tpg and cxj being the digital shearlets (see Definition 2.2) for j = ji,, . . ., ju, 
s £ {-2 j , • • • , 2 J }, and n€K,m€ U jtS . Then the FDST, which we abbreviate by S, 
takes the form 



5* = WJwP. 



(3.1) 

To be more precise, letting / be an image of size N x N and J w := y/wPI be the 
weighted pseudo-polar Fourier transform of /, the set of shearlet coefficients of / 
generated by the FDST can be written as 



SI = {eg : L ,n} U {ct- : L,j,s,m}, 



(3.2) 



where c L ° = (J w ,<pg)ci R for i = 0, 1 and c L j s m = (J w ,<Tj >s<m )Q R for i = 11, 12, 21, 22. 

Aiming to derive a closed form of S -1 , let P* denote the adjoint operator of P, 
which, for a function J : — > C, is given by 



P*J(u,v)= E J(u x ,u} y )e m ° 

(cj x ,w y )Gf2ij 



, U,V = -f, 



2 i- 
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Further, let W* denote the adjoint operator of W, which, given a sequence of shearlet 
coefficients Cj = SI as in (13. 21), is defined by 



Then we have the following result, which shows that the inverse FDST equals the 



adjoint FDST provided the weight function w satisfies (2.101 



Proposition 3.1. If w satisfies (2.101, then S 1 = P*^/wW* 



Proof. By hypothesis on w, we have P*wP — Id. Moreover, W*W = Id by 



Theorem 2.3 Since S = W^/wP by (3.1), the result is proved. □ 



Let us now assume that the weight function w does not fulfill the requirements 



in (2.10), i.e., P*wP ^ Id, which prevents simply using the adjoint operator. In this 
case, we mention two methods to compute the inverse. The first method is a direct 
approach by resampling trigonometric polynomials from the pseudo-polar grid to the 
Cartesian grid, which is discussed in detail in [Tj. The second method is an iterative 
approach by using the conjugate gradient method, which we now describe. 

Suppose we are given a sequence of shearlet coefficients Cj of some image /. By 



Theorem 2.3 the inverse digital shearlet windowing of / is y/wPI — W*Cj =: J, 



Consequently, the original image I can be computed by solving 



min II JwPI — Jui 2- 
Solving this problem is equivalent to solving the linear system of equations 

P*wPI = P*y/wJ w . (3.3) 
This shows that w plays the role of a preconditioner for the following normal equations: 

P*PI = P*J W . 



We refer to Subsection |2.2| for a selection of choices for weight functions and a dis- 
cussion about their performance as preconditioncrs. Since the matrix corresponding 
to P*wP is symmetric and positive definite, the conjugate gradient method can be 
used to solve (3.3). Algorithmic details will be discussed in Section [5] Let us just 



mention two main issues: The number of iterations required by the conjugate gradient 
method depends on the condition number of P*wP. Moreover, the conjugate gradient 
method only requires applications of P* and P to a vector, which can be computed 
0(N 2 log N) flops, in contrast to forming the complete matrices. 

4. Mathematical Properties of the FDST. Results on decay properties of 
the discrete shearlet coefficients and even more its sparse approximation properties are 
well-known, see [TS1 [231 US] . These continuum domain results do however not directly 
imply similar statements for the introduced digital shearlets due to the fundamentally 
different nature of a digital grid. Therefore, in this section, we will analyze decay 
properties of the digital shearlet coefficients for linear singularities. Moreover, we will 
prove shear invariance of the FDST. 

4.1. Decay Properties of FDST Coefficients for Linear Singularities. 

One main advantage of shearlets over wavelets is their ability to precisely resolve 
curvilinear - hence, in particular, linear - singularities due to their anisotropic shape 
and their additional direction-sensitive shear parameter |24j . Our computations will 
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show that this property carries over to the digital setting. Our model for a linear 
singularity will be a line. We remark that similar results can be shown for a Heaviside 
functions as a model. 

Let I be an image of size NxN with an edge through the origin of slope t satisfying 



|i| < 1, i.e., I(u,v) = 5(tu - v), -N/2 < u,v < N/2 - 1. Figure 4.1 indicates that 





Fig. 4.1. From left to right: A horizontal line in spatial domain and its pseudo-polar Fourier 
transform, and a line with slope < 1 in spatial domain and its pseudo-polar Fourier transform. 

for * 7^ an aliasing effect occurs. This effect is quite small in comparison with 
the intensity of the line, wherefore we will ignore it in the calculations. Also, by 
symmetry, we can assume * > without loss of gener ality. Setting m = |(EiV+ 1), 
for {u) x ,w y ) = (i£,-^f ■ §) € Q% we have 

N/2-1 N/2-1 
u=-N/2 u=-N/2 

If k = 0, we obtain I(uj X) U) y ) = N. If k ^ 0, 



«**> g-zit) fBn(irNk(l-%t)/(RN+i)) 



sin(7rfc(l- §t)/(RN + l)) 



Concluding, 



l(U) X ,U)y) = N, 



\I(w x ,U)y)\ 



< 



RN + 1 



2\k\\l 



t 



for k = 
for |fc| = 1,.. 



RN 
2 ■ 



Similarly, for (uj x ,uj v ) = (-§ • § , § ) G fijj, 



l(UJ X ,LOy) = N, 



\l(U X ,U)y)\ 



< 



(RN + l)N 



for I = t 



, for |fc| = 1, 



AT 



RN 



N 
2 ' 



We now distinguish two cases. On f^ 1 (similarly on fifj), for those digital shear- 
lets not on the seam lines, i.e., s £ {— 2 J ',2 J '}, we obtain 



21 \ I < 



E 



\I(u X ,Uy)\ 



J' s| (w^wjesupp tr^ o 

2- J '- 1 JV(s+l) 



< 



E 



iV/2 



4 J + 1 fl/2 

E 



2|fc| 
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2- : >- 1 N(s+t) 



E 



< 



< 



(RN + 1)N 

(RN + 1)N 
(RN + 1)N 1 



4 3 + 1 i?/2 

E 



I7V/2 - U\ 

l=2-i- 1 N(s-l) 1 ' 1 k=4i- 1 R/2 

2~ : >- 1 N(s+l) 



1 



E 



iV/2(l-2-J(s + l)i) 



4 3 + 1 i?./2 

E 

fe=4J'- 1 _R/2 



4^~[ 23- (1 -«) 
On Qjj 1 (similarly on fl R 2 ), 
1 



41og(2). 



1 

(4.1) 



1(^4 



s,m/ i'i? I 



< 



< 



1 



E 1-^(^)^)1 

(^ x ,w H )eSUpp crj^o 
2- J '- 1 A'(s+l) 4 j + 1 i?,/2 

E E (*(- 

t=2-i- 1 N(s-l) k=Ao- 1 R/2 



+ (l-5(s-t- f)5(£-t- N/2)) 



t- f)S(£-t- f)iV 
(RN + 1)N 



A\k\\e-t- f 



< 



(i&V + l)iV 



2" J " 1 Ar(s+l) 

E 



i 



4 .H 



r + 1 



E 



£=2-J- 1 Af(s-l),^^t-^ 



JV 1 71 fc=43-!i?./2 



5(a - i • f ) • (4-> +1 • f + 1)JV 



< (i?jV + l)jVlog(jV/2)41og(2) + <5( S -t-f)(4^+ 1 -f + l)A^ (i 



4^ 



j.s 



Digital shearlets on the seam lines can be dealt with similarly. Thus, by (4.1) and 
(4.2 1, we have upper bounds for shearlet coefficients which as j — > oo have asymptotic 
decay 

• 0(2- 3j / 2 ), if the shearlet is in and is not aligned with the line singularity, 

• 0(2~^ 2 ), if the shearlet is in 0]j and is not aligned with the line singularity, 

• 0(2 3 ^ 2 ), if the shearlet is in fl 1 ^ and is aligned with the line singularity. 

We now aim to show that line singularities can indeed be detected in the transform 
domain. For this, we will show a lower bound on the decay of the shearlet coefficients 
if the shearlet is aligned with the line singularity. To avoid unnecessary technicalities, 
we now assume I(u) x ,u> y ) = N for uj v /lj x = — 1/t and I(oj x ,u> y ) = for to y /tu x =/= — 1/t 
with < t < 1. The case for t < can be handled similarly. Notice that this condition 
merely assumes that the aliasing effects are negligible for the decay analysis. We 
further assume that the window function W and the bump function V are positive 
on their supports. On fijj, 

max{|(7,cr^ )n R |} > \{I, <r > )0 )n R | 

l 



J2 i(u x ,w y )W(A-iu x )Vj(s + 2i\ 
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' J2 S(s-t-f)-N-W(4-^ x )V^(s + 2^) 



^(s-f-f)-iV 

where Sj lS = | u] y )en R W(4~ : 'uj x )V(s + 2 J ^)| denotes the area of the window 
function with respect to j and s. Since Sj. s is approximately the same order as |7\Lj, s |, 
we conclude that 

N S 1 ■ 

maxiKl^wJa*]} > ~j^= = 0(V' 2 ) as j -> oo. 



This estimate combined with (4.2 1 shows that the asymptotic decay as j — ► oo in Q} R 
is 

• 0(2~ J / 2 ), if the shearlet is not aligned with the line singularity, 

• il(2 J / 2 ), if the shearlet is aligned with the line singularity. 

In this sense there is a strong difference between the decay rates of shearlet co- 
efficients between those aligned with the line singularity and those not aligned with 
the line singularity. 

Experimental results strongly support this analysis. Table [4~T] presents the max- 
imal absolute values of shearlet coefficients of an 512 x 512-image with a horizontal 
line (t = 0) and different scales, both aligned with the line singularity (c^^), and not 



aligned with the line singularity {c\, lax ). The data in Table 4.1 does clearly indicate a 
significant difference of decay rates between these two classes. We shall confirm this 
behavior from a different viewpoint, more precisely, a particular quantitative measure, 
in Subsection 17.51 

Table 4.1 

Maximal coefficients of a horizontal line 



3 


-l 





1 


2 


3 


4 


5 


c° 


0.160 


0.098 


0.068 


0.038 


0.018 


0.013 


2.6E-3 


c 1 


0.149 


0.048 


0.007 


1.9E-3 


4.8E-4 


2.3E-4 


2.4E-5 



4.2. Shear Invariance of the FDST. Since the shearing operator is quite 
distinctive in the definition of shearlets, one might ask whether the FDST is in fact 
even shear invariant. In the continuum setting, when choosing ip to be a shearlet 



generator as defined in ( 1.1 1, one can easily verify that, for any / e L 2 (R 2 ), 

(2 3 ^ 2 V;(^7 1 A 43 ■ -m),f(S s -)) = {2^ 2 iP{S~l 2] A A3 ■ -m),f). 

This identity can be viewed as a manifestation of shear invariance, since it states 
that the shearlet coefficient of a sheared image f(S s -) at scale j, shear k, and spatial 
position m equals the shearlet coefficient for the original image / at the same scale j, 
same spatial position m, but with a shear parameter shifted by 2- 7 s. 

The following results shows that the FDST associated with digital shearlets has 
a similar property. 

Proposition 4.1. Let I be an N x N image . Let j be a scale, s be a shear, 
t with \t\ < 1 be a slope, and m — (mi, m-i) be a spatial position such that 2 J i e Z 
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and —2 1 < s,s + 2H < 2 J . Let I t := I(Sf) be the sheared image of I such that 
i t (u) x ,u y ) = I(u> x ,Lo y - toj x ) for all (u x ,0J y ) G supp <T$ tBtm . Then 



(It, 0~j sm JQ R — (I, °~j.s+2i t.mr&R 



3 — 1 



Proof. It is sufficient to prove the claim for the case i = 21, since the other cones 
can be handled similarly. In this case, (u> x ,u> y ) = (^r,^f • ~^ L ) and hence 

f~ ~. \ _ l + \ _ l2k 2k -2{l+tN/2)\ 

[L0 X , U)y) :— [U X ,LOy — UO X ) — {^ft, -ft ■ jy )■ 

Since I t = I((S^ ) T -) on the support of s m , we have 



Tlj,s\ ■ (h,o- j s m )n R 
= J2 it(u x ,u y )W(4-iLJ x )Vi( S + 27|^) e - 2 -(™.( 4 " 3 -- 23 '^)) 

5^ 7(^,w y )I^(4-^ :c )^(.s + 2^ + 2^)e- 2 "( m ^ 4 "^' 23 S)) 
(s a ,,o y )e(s t - 1 )i-ng 

= fa". 0y)W(4-i«.) + 2*t + a, gLje-^*-'^ fe+ a '*)) 

(a> e ,a 1) )e(sr 1 )' r "li 



This proves the claim. □ 

5. Implementation of the FDST and its Inverse. After the formal intro- 
duction of the FDST and its theoretical analysis, we now turn to discuss the details 
of our implementation, including the forward transform, the adjoint transform, and 
the inverse transform. The associated code ShearLab-PPFT-1 . can be downloaded 
from www . ShearLabTorg[ 



Figure 5.1 provides an overview of the main steps of of the FDST and its inverse, 
which were defined in Sections [2] and [3j respectively. 



Image I 








Peudo-Pokr Image J 








Weighted Image J» 








Shear|el Coefficients: 



Image I 








Peudc-Polar linage J 








Weighted Image J* 








Shcar|et Coefficient: 
Cj.k.m 



Adjoint PPFTor 
Inverse Using CG 



Fig. 5.1. Flowcharts of the FDST (left) and its inverse (right). 
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5.1. Choices of Parameters. Before discussing the implementation details of 
each transform, let us start by elaborating on the choices of the main parameters 
involved in our design of the FDST. 

5.1.1. Choice of the Parameter tuq. The definition of the pseudo-polar Fou- 
rier transform (2.5| has Too as a free parameter. In the paper pQ, in which only 
the oversampling rate R = 2 was considered, mo was chosen to be 2N + 1. This 
particular choice enabled utilization of the 1D-FFT for the implementation of the 
fast pseudo-polar Fourier transform. In the situation of an arbitrary oversampling 
rate R, which we consider, the fast pseudo-polar Fourier transform is based on a 
(discrete) fractional Fourier transform. In order to utilize the 1D-FFT along one 
direction, m n has necessarily to be chosen as j^(RN +1); for details see Subsection 
|5.2| When choosing a different parameter too, the fractional Fourier transform needs 
to be applied on both directions, which certainly would significantly lower the speed 
of the FDST. The (discrete) fractional Fourier transform can be implemented with 
the same complexity as the 1D-FFT, but a different constant; in fact it is about 5 
times slower than the 1D-FFT, see [2]. To accelerate the speed, in our ShearLab 
package, the default too is consequently set to be ji(RN + 1). 

5.1.2. Choice of Weight Function w. A second free parameter is the weight 
function w, for which the criterion for isometry in Theorem |2.1| is required; and we 
discussed some choices in Subsection |2.2.2| As could be seen, the performance in 
terms of almost isometry Mi Som differs depending on the types of images which are 
considered. Hence the weight needs to chosen depending on the application. It should 
be also emphasized that the particular type of weighting we considered in Subsection 
|2.2.2| provide good preconditioners for the conjugate gradient method, in case an even 
more accurate inverse than the adjoint is required. In our ShearLab package, various 
choices of w are available. 

5.1.3. Choice of Window Functions Wo,Vo and W,V. The final main pa- 
rameter is the choice of the window functions Wo, Vq and W, V. In our implementation, 
we use Meyer wavelets and define Wo and W to be the Fourier transform of the Meyer 
scaling function and wavelet function, respectively, i.e., 



1 

cos[fi/(f|£| 




§)] 



iei< 



i 

4 ' 



!<l£l<i, 

otherwise, 



and 



w(0 



< 



< 1, 



i < ici < 4, 

otherwise, 



where v is a C function or C°° function such that 



v{x) 



v{\ - x) 



x < 0, 
< x < 1, 
x > 1. 



(5.1) 



In ShearLab, v is set to be v(x) — 2x 2 for < x < 1/2 and v{x) = 1 — 2(1 — x) 2 for 
1/2 < x < 1, which is a C 1 function. Other choices are also available. The choice of v 
then fixes Wq and W. Since |Wq(£)| 2 + |W(£)| 2 = 1 for |£| < 1, the required condition 



(2.13) is satisfied. These choices for Wo, W, and v are illustrated in Figure 5.2 
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Provided v satisfies (5.1 ), it can be used to design the bump function V by setting 



= a/z^I +£) + v{l - t) for -1 < f < 1. This automatically satisfies ( |2.15[ ), and 
is our choice in the implementation. The function Vo is defined to be Vq = 1. 

5.2. FDST. We now present our implementation of the FDST, which mainly 
consists of three parts: The fast pseudo-polar Fourier transform, the weighting on 
the pseudo-polar grid, and the digital shearlet windowing on the pseudo-polar grid 
followed by 2D-1FFT. These parts correspond to the operators P, w, and W intro- 
duced and discussed in Section[3] To facilitate our presentation of the implementation 
details, we require the following operators: 

The (unaliased) fractional Fourier transform (frFT) of a vector c G C JV+1 with 
respect to a fraction a € C is defined to be 



N/2 

(F N+1 c)(k):= £ c(j)e- 2 ^ k - a , fc = -f,. 

j = -N/2 



N 
2 ■ 



It was shown in [2], that the fractional Fourier transform F^ +1 c can be computed 
using 0(N log N) operations. In the special case of a = 1/(N + 1), the fractional 
Fourier transform becomes the (unaliased) ID Fast Fourier Transform (1D-FFT), 
which in the sequel we will denote by F±. Similarly, we will denote the 2D Fast 
Fourier Transform (2D-FFT) by F 2 , and the inverse of the F 2 by F 2 _1 (2D-1FFT). 

Let now N be even and m > N be an odd integer. Then the padding operator 
Em,n acting on a vector c € C N gives a symmetrically zero padding version of c in 
the sense that 



{E m<N c)(k) 




K _ 1 
' 2 ' 

m\\r N_ N i\ 

• ) 2 J \ l 2 ' ' " ' ' 2 /' 



The fast pseudo-polar Fourier transform was introduced in pQ, however only for 
the oversampling rate R = 2. We will next extend this algorithm to an arbitrary 
oversampling rate and show that its complexity is also 0(N 2 log N). For this, let I 
be an image of size N x N. We restrict to the cone Q R . Choosing mo = ^(RN + 1) 
to utilize the 1D-FFT, for (u x ,u y ) € Q^, 

N/2-1 N/2-1 N/2-1 

/(w bjWj ,)= ]T I(u,v)e-^ {u ^ +v ^ ) - E I(u,v)e- 2 ^ {u ^ +v ^ ) 

u,v=-N/2 u=-N/2v=-N/2 
N/2-1 

x I(u,v)e e (r«+i).iv. 
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The above identity shows that / on il]^ can be obtained by performing 1D-FFT on 
the extension of I along direction v and then applying the frFT along direction u. 
More precisely, let I(u, ■) := Ern+i,nI(u, ■), —N/2 < u < N/2 — 1 be the symmetric 
zero padding of / and let I\(u, •) := Fil(u,-) be the ID FFT of I along direction 
v. Also, let I%(-,k) := En+i,n^i{-, k) for k = — ^,...,4^ be the symmetric zero 
padding of 1\ along direction u. Note that, 1\ is then of size (RN + 1) x N and I 2 is 
of size (RN + 1) x (N + l).Then the above identity can be written as 

N/2-1 N/2 

I(uj x ,u) y )= 7 1 (ii,fc)e~ 27r ™^(Wi' C )-«/2 = l2(u,k)e- 2M (™+i>" 

u=-N/2 u=-N/2 

= (F^ +1 U,k))(i), 

where ctk = — ^hm+i)n/2 * s ^ e fraction in the frFT. 

Thus, the pseudo-polar Fourier transform I(u) x ,cj y ), (oj x ,ujy) £ Q, 1 ^ - similarly for 
fi^j - for arbitrary oversampling rate R can be computed by the steps described in 
Algorithm [T] 



Algorithm 1 Fast Pseudo-Polar Fourier Transform of / on fl^ 

(a) Input: Image I of size N x N. 

(b) Output: The pseudo-polar Fourier transform I^i . 

(c) PPFT: Image I(u, v),-N/2 < u, v < N/2 - 1. 

1: I(u,-) Ern+i,nI(u, ), —A/2 < u < N/2 — 1. Symmetrically padding the 
image I along direction v to obtain an image / of size (RN + 1) x A. 

2: I\(u, ■) F\I(u, •), —N/2 < u < N/2 — 1. For each vector of I along direction v, 
perform the 1D-FFT along direction v to get I\. 

3: I 2 (-,k) <- E N+hN h(-,k), -RN/2 < k < RN/2. Symmetrically padding the 
image 1\ along direction u to get image I2 whose size along direction u is N + 1. 

4: I q i r <- F^" +1 I 2 (-, k), -RN/2 <k< RN/2. Perform the fractional Fourier trans- 
form with respect to a fraction — — rgj^upfn a l° n S direction u. 



Since the 1D-FFT and lD-frFT require only 0(N log A) operations for a vector 
of size A, the total complexity of the pseudo-polar Fourier transform is 0(N 2 log A) 
for an image of size A x A. It should be emphasized that our Algorithm 1 is more 
simple and efficient even in case R = 2 than the algorithm described in pQ, in which 
F 2 is applied to / followed by the application of F-f 1 to the resulting image in order to 
obtain I\ . The algorithm in [1 is slightly more redundant while ours is more efficient 
by combining these two steps to one step. Moreover, the fractional Fourier transform 
along direction u is only of size A + 1, as compared to size RN + 1 in pQ. 

We would like to also remark that, for a different choice of constant mo, one 
can compute the pseudo-polar Fourier transform also with complexity 0(A 2 logA) 
for an image of size A x A, in which case the fractional Fourier transform needs to 
be applied in both directions u and v of the image. For example, for computing I 
on f2^, the fractional Fourier transform is first applied to I along direction v with a 
fixed fraction constant a — — ^tttt , and then to the another direction u with fractions 

moH/2 ' 

afc = mn fl/2-iv/2 depending on k. In this case, it is, of course, slower than the special 
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choice mo = jj{RN + I), since frFTs are involved in both directions while the special 
choice mo utilizes FFT. 

The next step, which is an application of the weight function, is simply a point- 
wise multiplication on the pseudo-polar grid, i.e., 

(d) Weighting: Let J : fi/j — > C be the pseudo-polar image of I and w : 
CIr — > K + be any suitable weight function on flji. Compute the point-wise 
multiplication J w (u x ,oj y ) = J(uj x , u> y ) ■ y/w(oj x ,ujy), (to x ,uJy) € CIr. 

For the windowing, a sequence of subbands {</? ° : L o}^Wj s o : h s i L } i s computed, 
whose frequency tiling covers the pseudo-polar grid. As discussed in Subsections |2.3| 
and |5.1| these subbands are constructed using Meyer scaling functions and wavelets. 
The weighted image is then windowed by each subband, followed by an application of 
the 2D-iFFT to each windowed subimage. This results in a sequence of coefficients - 
the digital shearlet coefficients - structured block by block: 

(e) Subband Windowing: Compute the digital shearlet coefficients by window- 
ing J w with respect to the digital shearlets {</?^° : l , n} U {crlj m : j, s, m, t} 
defined in Definition 12.21 

For the detailed step-by-step description of the algorithm for the FDST used in 



ShearLab 



we refer to Algorithm A.f (see Appendix) 



5.3. Adjoint FDST. The algorithm for the adjoint shearlet transform can be 
straightforwardly derived from the FDST, and the main idea was already described 
in Section [3] It should be noted that the adjoint fractional Fourier transform for a 
vector c G C N+1 with respect to a constant a € C is given by F^^c. Moreover, 
for m > N the adjoint operator E* n N for the padding operator E m .N is given by 
(E* mN c)(k) = c(fc), k = -N/2, . . . , N/2 - 1 for a vector c € C m . The adjoint shearlet 
transform can be computed with a complexity 0(N 2 log N) similar to the FDST, since 
it is obtained simply by 'running the FDST backwards'. 

For the detailed step-by-step description of the algorithm for the adjoint FDST 
used in ShearLab, we refer to Algorithm |A.2| 

5.4. Inverse FDST. The main idea to use the conjugate gradient method was 
already described in Section [3j and for a detailed step-by-step description of the 



algorithm for the inverse FDST used in ShearLab, we refer to Algorithm A. 3 



6. Quality Measures for Algorithmic Realization. To ensure and also 
prove that our implementation satisfies the previously proposed desiderata, we will 
now define quality measures for each of those and in the sequel provide numerical 
results on how accurate our implementation satisfies these. It is moreover our hope 
that these measures shall serve as comparison measures for future implementations. 
Although some measures are stated in 'shearlet language', most are applicable to any 
directional transform based on parabolic scaling. 

In the following, P shall denote the pseudo-polar Fourier transform defined by 
Algorithm [T]), w shall denote the weighting applied to the values on the pseudo-polar 
grid, W shall be the windowing with additional 2D-iFFT, S shall denote the FDST 
defined by Algorithm |A.1[ and S* its adjoint defined by Algorithm A. 2 We will 



further use the notation GaJ for the solution of a matrix problem AI = J using 
conjugate gradient method with residual error set to be 10 -6 . 

[Dl] Algebraic Exactness. 

Comments: We require the transform to be the precise implementation of a theory 
for digital data on a pseudo-polar grid. In addition, to ensure numerical accuracy, we 
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provide the following test, which provides a quantitative measure for the closeness of 
the windows to form a tight frame. 

Measure: Generate a sequence of 5 random images Ii , . . . , 1$ (the integer 5 is chosen 
for the purpose of fixing the number of images to enable precise comparison) on the 
pseudo-polar grid for N = 512 and R = 8 with standard normally distributed entries. 
Our quality measure will then be the Monte Carlo estimate for the operator norm 
||W*W-M||op given by 

w \\W*Wh-h\\ 2 
M aig = max — . 

»=1,...,5 \\Ii\\ 2 
[D2] ISOMETRY OF PSEUDO-POLAR FOURIER TRANSFORM. 

Comments: To ensure isometry of the utilized pseudo-polar Fourier transform, we 
introduced a careful weighting of the pseudo-polar grid. The following test will now 
measure the closeness to being an isometry. We expect to see a trade-off between 
the oversampling rate and the closeness to being an isometry. Since the measure 
shall however serve as a common ground to compare different algorithms based on 
parabolic scaling, we do not take the oversampling rate into account in the proposed 
measure. Instead, we would like to remind the reader that this rate will instead 
affect the measure for speed [D6]. In the sequel, we will now provide three different 
measures, each being designed to test a different aspect. 

Measure: 

• Closeness to tightness. Generate a sequence of 5 random images I\ , . . . , I§ of 
size 512 x 512 with standard normally distributed entries. Our quality measure 
will then be the Monte Carlo estimate for the operator norm \\P*wP — Id|| op 
given by 

M - max W P ^ P ^-^h 

lvl isomi — llldA 11 7- II 

»=1,...,5 ||Vi|| 2 

• Quality of preconditioning. Our quality measure will be the spread of the 
eigenvalues of the Gram operator P*wP given by 

M _KUpyp) 

M isom2 - Xm . n{p * wp) ■ 

• Invertibility. Our quality measure will be the Monte Carlo estimate for the 
invertibility of the operator yfwP using conjugate gradient method G^ P 
given by 

WG^py/wPIj - Ij\\ 2 

M isom 3 — . max c u T u 
»=1,...,5 \\Ii\\ 2 

[D3] Tight Frame Property. 

Comments: We now combine [Dl] and [D2] to allow comparison with other transforms, 
since those singleton tests might not be possible for any transform due to a different 
inner structure. 

Measure: Generate a sequence of 5 random images I\, . . . ,1$ of size 512 x 512 with 
standard normally distributed entries. 
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• Adjoint transform. Our quality measure will be the Monte Carlo estimate for 
the operator norm \\S*S — ld\\ op given by 

\\S*SIi- Ii\\ 2 
M tightl = max — jj . 

t=l,...,5 ||7j|| 2 

• Inverse transform. Our quality measure will be the Monte Carlo estimate for 
the operator norm \\G^ P W* S — Id|| op given by 

\\G^ P W*Sh-h\\ 2 
M tig ht 2 = max, iTFli ' 

»=1,.. .,5 llijl^ 

[D4] Space-Frequency-Localization. 

Comments: The purpose of this test is to provide quantitative measures for the degree 
to which the windows - in our case the digital shearlets - are localized in both space 
and frequency. For this, we numerically measure mathematically precise notions of 
both decay and smoothness. 

Measure: Let I be a shearlet in a 512 x 512 image centered at the origin (257, 257) 
with slope s = at scale j — 3, which in our case is the shearlet 03 0,0 + a 3fl,o- Our 
quality measure will then be four-fold: 

• Decay in spatial domain. Compute the decay rates d\, . . . ,d 5 i 2 along lines 
parallel to the y-axis starting from the line [257, : ] and the decay rates 
(^512, • • ., c?io24 with x and y interchanged in the following way: Consider 
exemplarily the line [257 : 512,1]. First compute the smallest monotone 
majorant M(x, 1), x — 257, . . . , 512 - note that we could have also chosen a 
different 'envelope' - for the curve \I(x, 1)|, x = 257, . . . , 512. Then the decay 
rate along the line [257 : 512, 1] is defined to be the average slope of the line, 
which is a least square fit to the curve log(M(x, 1)), x — 257, . . . , 512. Based 
on these decay rates, we choose our measure to be the average of the decay 
rates given by 

M d eca yi = ^7^7 ^ ' 

i=l,. ..,1024 

• Decay in frequency domain. To check whether the Fourier transform of / is 
compactly supported and to check its decay rate, let I be the 2D-FFT of I 
and compute the decay rates di, i = 1, . . . , 1024 as before. Then we define 
the following two measures: 

O Compactly supportedness. 

M _ max H , H < 3 \i(u,v)\ 
max Hi „ \I{u,v)\ 

O Decay rate. 

M decay2 = — — - di. 

i=l,...,512 

• Smoothness in spatial domain. Smoothness will be measured by the average 
of local Holder regularity. For this, for each (u , vq) e {1, . . . , 512} 2 , compute 
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M(u, v) = \I(u, v) — I(uq, i>q)|, < max{|M — Uo\, \v — Vq\} < 4. Then the 
local Holder regularity ot UOtVo is the least square fit to the curve (u,v) t— > 
log(\M(u,v)\). The measure is given by 

Msmoothi — -^o? ^ ; C^u.v • 

• Smoothness in frequency domain. Smoothness will again be measured by the 
average of local Holder regularity. Proceed as for measuring the smoothness 
in spatial domain now applied to j, the 2D-FFT of I, to derive the local 
Holder regularity a U) „ for each (u,v) G {1, . . . ,512} 2 . The measure is then 
given by 

u,v 

[D5] True Shear Invariance. 

Comments: This test shall provide a measure for how close the transform is to being 
shear invariant. In our case, the theory gives 

(2^/ 2 i,(S^A 4] ■ -m), /(&•)) = (2 3 ^ 2 i,(S-l 2]s A 4J ■ -m),f) , 

and we expect to see this behavior in the shearlet coefficient as discussed in Subsec- 
tion [O 



Measure: Let / be an 256 x 256 image with an edge through the origin (129, 129) of 
slope 0. Fix s = 1/2, generate an image I s :— I(S S -), and let j £ {1,2,3,4} be a 
scale. Note that 2 J s G Z for each of j € {1, 2, 3, 4}. Our quality measure will then be 
the curve 

M shear j = max — f , scale j = 1, 2, 3, 4, 

-23<fc,fc+23s<2J 1 1 7 1 1 2 

where Cj k is the set of coefficients - in our case the shearlet coefficient - at scale j 
and at shear index k. 

[D6] Speed. 

Comments: When testing the speed, not only the asymptotic behavior will be mea- 
sured but also the involved constants. It should further be mentioned that for practical 
purposes it is usually sufficient to analyze the speed up to a size of N = 512. 

Measure: Generate a sequence of 5 random images 7,, i = 5, . . . , 9 of size 2 l x 2* with 
standard normally distributed entries. Let Si be the speed of the transform S applied 
to Ii. Our hypothesis is that the speed behaves like s; — c ■ (2 2l ) d ; 2 21 being the size 
of the input. Let now d a be the average slope of the line, which is a least square fit 
to the curve i i-> log(si). Let also be the 2D-FFT applied to Ii, i = 5, .. .,9. Our 
quality measure will then be three-fold: 
• Complexity. 
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The constant. 



1^ 



Comparison with 2D-FFT. 



5 U ( 221 ) 



i 9 

i=5 ^ 

[D7] Geometric Exactness. 

Comments: Geometric objects such as edges should be as precise as possible be resem- 
bled by the coefficients in the sense that analyzing the decay of the coefficient should 
detect such features. For the FDST, these properties were theoretically analyzed in 
Subsection 4.1 Our model will be an image containing one line of a particular slope. 

Measure: Let be 256 x 256 images of an edge through the origin (129, 129) 

and of slope [— 1, —0.5, 0, 0.5, 1] and the transpose of the middle three, and let c,j- be 
the associated shearlet coefficients for image at scale j. Our quality measure will 
be two-fold: 

• Decay of significant coefficients. Consider the curve 

8 



1 



1 = 1 



max \aj(oi shearlets aligned with the line) | , scale j, 



let d be the average slope of the line, which is a least square fit to the logarithm 
of this curve, and define 

Mg eoi — d. 

• Decay of insignificant coefficients. Consider the curve 

1 8 

- max |cj j (of all other shearlets) | , scale j, 
8 r— f 

let d be the average slope of the line, which is a least square fit to the logarithm 
of this curve, and define 

M ge02 = d. 

[D8] Robustness. 

Comments: Two different types of robustness will be analyzed which we believe are 
the most common impacts on a sequence of transform coefficients. We wish to men- 
tion that we certainly also could have considered additional manipulations of the 
coefficients such as deletions; but to provide sufficiently many tests balanced with a 
reasonable testing time, we decided to restrict to those two. 

Measure: 

• Thresholding. Let I be the regular sampling of a Gaussian function with 
mean and variance 256 on {—128, 127} 2 generating an 256 x 256-image. 
The quality measure for k = 1, 2 will be the curve 



M, 



thres kiPk 



\Gy^pW* thres fc , Pfc SI - I\\ 2 
IUII 2 



where 
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— thresi jPl discards 100 • (1 — 2 Pl ) percent of the coefficients with pi — 
[2:2:10]), 

— thres2 lP2 sets all those coefficients to zero with absolute values below 
the threshold m/2 P2 with m being the maximal absolute value of all 
coefficients with p 2 = [0.001 : 0.01 : 0.05]). 

• Quantization. Let / be the regular sampling of a Gaussian function with 
mean and variance 256 on {—128, 127} 2 generating an 256 x 256-image. 
The quality measure will be the curve 

M WGy^pW* quant, SI - I\\ 2 _ _ 

M quant , q = — r — , q = [8 : -0.5 : 6], 

W 1 \\2 

where quant ? (c) = round(c/(m/2 9 )) • (m/2 9 ) and m being the maximal ab- 
solute value of all coefficients. 

7. Numerical Evidence. In this section, we provide numerical results for the 
tests [D1]-[D8] detailed in Section [6] of our present implementation, which we con- 
sider as having reached a mature state after careful tuning the parameters depending 
on these performance measures. The associated code ShearLab-PPFT-1 .0 can be 
downloaded from www . ShearLab . org| 



7.1. Results for Tests [Dl]— [D3]. Table 7.1 presents the performance with 
respect to the quantitative measures in [D1]-[D3]. 

Table 7.1 
Results for [D1]-[D3] 



M alg 










M t ight 2 


6.6E-16 


9.3E-4 


1.834 


3.3E-7 


9.9E-4 


3.8E-7 



The quantity M a i g sa 6.6E-16 confirms that the VSH defined in Definition 2.2 is 
indeed up to machine precision a tight frame. 

The slight tightness deficiency of M tightl sa 9.9E-4 (also M isomi s» 9.3E-4) mainly 
results from the isometry deficiency of the weighting. However, for practical purposes 
this transform can be still considered to be an isometry allowing the utilization of the 
adjoint as inverse transform. Progress on the choice of weights will further improve 
this measure. Observe though that there is a trade-off between the sophistication of 
the weights, the running time of S, and the smoothness of the shearlets. 

Further, note that the condition number (Mi Som2 sa 1.834) of the Gram matrix 
is quite close to 1, which - in case an even higher accurate inverse than the adjoint 
is required - allows us to employ the conjugate gradient method very efficiently for 
computing the inverse of the FDST (M isom3 sa 3.3E-7 and M Ught , 2 w 3.8E-7). 

7.2. Space- Frequency-Localization Test [D4]. The reference shearlet / re- 



quired for [D4] is illustrated in Figure 7.1 in both the spatial and frequency domain. 

Table 7.2 presents the space-frequency-localization measures Md ecayi , Mdecay 2 i 
and M supp for analyzing the decay in spatial and frequency domains as well as the 
measures M srnoo t/ l:1 and M smoot h 2 for smoothness in spatial and frequency domains. 

The measurements indicate that the shearlet illustrated in Figure |7.1| decays 
slower in spatial domain (Mj jecayi -1.920) than in the frequency domain (Md ecay2 sa 
-3.257) in terms of the average decay rates. In terms of average local Holder smooth- 
ness, the shearlet in spatial domain is smoother than in the frequency domain proven 
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Fig. 7.1. Graphs of the reference shearlet for test [D4] centered at origin at scale 3 and shear 
parameter in both spatial and frequency domain with N = 512. 

Table 7.2 
Results for [D4] 



decay i 




M decay 2 


^-smoothi 


M smooth,2 


-1.920 


5.5E-5 


-3.257 


1.319 


0.734 



by the fact that M smoothl w 1.319 > 0.734 w M smoottl2 . The very small value M supp 
5.5E-5 - although it is not zero due to round off errors - indicates that our shearlet 
is indeed compactly supported. 

7.3. Shear Invariance Test [D5]. Test [D5] requires a reference image / and 
a sheared version I s = I(S S -), which is illustrated in Figure 7.2 for s — 0.5 alongside 
with their pseudo-polar Fourier transforms I and I s , respectively. 




Fig. 7.2. Graphs of the reference image I and its sheared version I s for s = 0.5 in spatial 
domain and I and I s in pseudo-polar domain. 



The shear invariance measurements are presented in Table [773) This table shows 
that the FDST is indeed almost shear invariant. A closer inspection shows that 
M s hear,i and M s h ea r.2 are relatively small compared to the measurements with respect 
to finer scales M s h ea r,3 and M s hear,4- The reason for this is the aliasing effect which 
shifts some energy to the high frequency part near the bou ndar y away from the edge 
in the frequency domain, see also the graph of I s in Figure 



7.2 



7.4. Speed Test [D6]. The measurements for speed performance, i.e., for the 
complexity M speei i x , the constant M spee( i 2 , and the comparison with 2D-FFT M speec t 3 , 
are presented in Table |7.4| 

The complexity constant Mspe^ ~ 1-156 confirms that the order of complexity 
of our FDST algorithm is indeed close to linear (compare 0(N 2 log N) in theory). 
The constant of comparison with 2D-FFT is M spee d 3 ~ 280, which seems signifi- 
cantly slower than the usually 2D-FFT. However, we should notice that the 2D-FFT 
is applied directly to the image of size N x N while the FDST employs fractional 
Fourier transforms and subband windowing on a oversampling pseudo-polar grid of 
size 2 x (RN + 1) x (N + 1) with oversampling rate R = 8. Taking into account 
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Table 7.3 
Results for [D5] 



Scale 


l 


2 


3 


4 


M shear, j 


1.6E-5 


1.8E-4 


0.002 


0.003 






Table 7.4 










Results for [D6]-[D7] 








M speedy 


M speedy 


M ge0l 


Mg eQ2 


1.156 


9.3E-6 


280.560 


-1.358 


-2.032 



that the fractional Fourier transform is about 5 times slower than the FFT and that 
the redundancy (about 4i?) comes from subband windowing with oversampling rate 
R = 8, we conclude that the constant M speec i 3 is reasonable and our implementation 
is indeed comparable with the 2D-FFT. 

7.5. Geometric Exactness Test [D7]. Test [D7] requires reference images 
with edges at various slopes which are plotted in Figure [73] 




Fig. 7.3. The graphs of edges in both spatial and frequency domains, 
domain. Bottom 8: Edges in frequency domain. 



Top 8: Edges in spatial 



The associated geometric exactness measurements M ge0l and M ge02 are presented 



in Table 7.4 These two measurements show the differences between the decay rates 
of the shearlet coefficients aligned with edges (significant shearlet coefficients) and the 
shearlet coefficients not aligned with the edges (insignificant shearlet coefficients). As 
theoretically proven in Subsection |4.1| and hence expected for the implementation, 
the insignificant shearlet coefficients (M ge02 ps -2.032 ) decay much faster than the 
significant shearlet coefficients with increasing scale j. 

[D8]. Table 7.6 presents the measurements for the ro- 



7.6. Robustness Test 

bustness test [D8]. 

These results show that even if 100(1 — 2~ 10 ) pa 99.9% of the shearlet coefficients 
are discarded, the original image is still well approximated by the reconstructed image 
(Mthres! p ~ 0.007). Thus the number of the significant coefficients is relatively small 
compared to the total number of shearlet coefficients. The second row indicates that 
knowledge of the shearlet coefficients with absolute value greater than m(l — 1/2 001 ) - 
hence about 0.1% of coefficients - is sufficient for precise reconstruction (M t f lresi p2 
0.005). The quantization test M quant , q attests the FDST high resilience against even 
quite coarse quantization. 



SHEARLAB: A RATIONAL DESIGN OF A PARABOLIC SCALING ALGORITHM 35 

Table 7.5 
Results for [D8] 



M t hres ljPl 


1.5E-8 


7.2E-8 


2.5E-5 


0.001 


0.007 




0.005 


0.039 


0.078 


0.113 


0.154 


■Equant ,g 


0.034 


0.047 


0.057 


0.071 


0.109 



8. Conclusions. We described a natural digitization of the continuous shearlet 
transform based on the pseudo-polar Fourier transform. For this, we first introduced a 
digital shearlet tight frame, which is an exact digitization of the classical band-limited 
system defined in the continuum domain. This shows that the shearlet framework pro- 
vides a unified treatment of the continuum and digital realm. We then described a fast 
digital shearlet transform (FDST) whose main ingredients are a weighted pseudo-polar 
Fourier transform to achieve an isometric mapping into the pseudo-polar domain and 
a windowing based on the digital shearlet tight frame. We discussed several choices of 
weighting functions to achieve almost isometry so that inverse digital shearlet trans- 
form can be obtained by using the adjoint transform. Various properties of the FDST 
are proven, and details of its implementation are provided. We further define several 
measures that quantify the performance of the FDST, thereby justifying the term 
'rational design'. These measures also provide a common ground to compare the per- 
formance of different algorithms based on parabolic scaling. The FDST as well as the 
performance measures are publicly available in the software package ShearLab, see 
www . ShearLa b . org| 

Appendix A. Algorithms. 

In this appendix, we provide pseudo-codes in Matlab style for algorithms of FDST 



(Algorit hm [AJJ ), of the adjoint FDST (Algorithm |A2j), and the inverse FDST (Al- 
gorithm A.3|. For the required notation, we refer to Section [HI 
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Algorithm A.l Fast Digital Shearlet Transform (FDST) 

(a) Input: Image I — {[I] u ,v '■ —N/2 < u, v < N/2 — 1} of size N x N, the oversam- 
pling rate R, and the precomputed weight matrix w of size 2 x (RN + 1) x (N+ 1). 

(b) Output: Digital shearlet coefficients Cj = {&° : lq} U {cf- s : j, s, i}. Each s is 



a matrix of size £] x £| s , see Definition 



2.2 



(c) PPFT: The pseudo-polar Fourier transform (see also Algorithm [lj . 
1: J <- zeros(2, RN + 1, N+ 1). m <- RAT + 1. //J is the grid J = J n i U J Q 2 r . 

for sector = 1 to 2 do 

if sector = 2 then I = I T end if 
for v = -N/2 to iV/2 1 do 

9 i 9 ErN+1,NQ , [J]sector,:,v <~ Fiq. 

//symmetric zero-padding and FFT along direction u. 
end for 

for k = -RN/2 to RN/2 do 

q <— [</]sector,fc,:j ^ mN/2 ' 9 *~ ^N+l?> [^]sector,fc,: *~ Q- 

//frFT along direction u . 
9: end for 
10: end for 

(d) Weighting: Weighting on the pseudo-polar grid with precomputed weight w. 
11: J-^J.*y<w. //.* is the point -wise multiplication. 

(c) Subband Windowing: Subband windowing with T>SH. 
12: L <- -riog 4 (-R/4)"|, H <- riog 4 7V], sMax <- 2 ■ 2 H + 1. 
13: Ci <- cell{4, H-L + 2, sMax}. 
14: for l = 11,12,21,22 do 
15: for j = L to H do 

16: if j < then tile <- else tile 2 3 end if 
17: for s = —tile to tile do 

18: Let Cj s be a submatrix of J of size £] x £ ; 2 S with respect to the support 

of the digital shearlet a 1 - . 

Cj, s <" 4«-*°U^> c j, s <" ^T^,.- //Windowing with 2D iFFT. 

C/{t,i,s} <- & j s . 

end for 
end for 
end for 

Let ip\, ifQ be the shear lets associated with the low- frequency part. Let c 4 , i = 1,2 
denote the submatrix of J with respect to the support of </?q, i = 1,2. 

25: C* «- c\*<^, Z = 1,2. 

26: Cj{i,L- 1,0} 4- c l , i = 1,2. 
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Algorithm A. 2 Fast Adjoint Digital Shearlet Transform (Adjoint FDST) 

(a) Input: Digital shearlet coefficients Ci = {c L ° : l} U {cj : j, s, i}, where each cj 
is a matrix of size £] x s (see Definition 2.2 1, the oversampling rate R, and the 
precomputed weight matrix w of size 2 x (7Z7v + 1) x (A + 1). 

(b) Output: An image / = {[!]„,„ : -A/2 < u,u < A/2 - 1} of size A x A. 

(c) Adjoint Subband Windowing: Subband windowing by VSH. 

h L < |"log 4 (i£/4)~|, <- riog 4 7V], J«-zeros(2,.RA + l,A + l). 

2: c* <- C/{i,L- 1,0}, z = 1,2. 

3: C l <r- C l .*(fi l Q , 1 = 1,2. 

4: Assign c 1 to J with respect to the support of ip},, i = 1,2. 
5: for i = 11,12,21,22 do 
for j = L to H do 

if j < then tile <— else tiZe <— 2 3 end if 
for s = —tile to iiZe do 



c: 



3,s 



<- Ci{lJ, s}. 

<- F 2 c' j s ./ /the 2D-FFT. 



c j,s *~ c j,s-* a j,s,0- 

Assign Cj. s to J with respect to the support of crj s0 . 
end for 
end for 
end for 



(d) Weighting: Weighting with the precomputed weight w. 
16: J-^J.*^/w. //■* is the point -wise multiplication. 

(c) Adjoint PPFT: The adjoint pseudo-polar Fourier transform. 
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I <- zeros{N, A), I Q <- zeros(N, N),m <- RN + 1, J x «- zeros(m, A), 
for sector = 1 to 2 do 

for k = -RN/2 to RN/2 do 

9 ^ [J] sector, k,: ■ & mN/2 ' 9 ^ ^N+lQ' 

«- E^ +1N q. 

end for 

for w = -AT/2 to A/2 - 1 do 

q <- [Ji]-.,v q <- Fiq, q <- E* n N q. 

[h]:,v <~ ?■ 

end for 

I = I + I . 
end for 
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Algorithm A. 3 Inverse Fast Digital Shearlet Transform (Inverse FDST) 

(a) Input: Digital shearlet coefficients Cj = {c L ° : to} U {c) s ■ j,s,t}, the oversam- 
pling rate R, the precomputed weight matrix w of size 2 x (RN + I) x (N + 1), an 
initial guess Iq, a precision parameter s, and a maximal iteration number it Max. 

(b) Output: An image I = {[I] u , v : -N/2 < u, v < N/2 - 1} of size N x N. 

(c) Initialization: Perform the adjoint shearlet windowing on Ci by using Algo- 
rithm 2 to obtain an image 6, say. 

(d) CG Iteration: Let A = P*wP generated by Algorithm]!] Algorithm A.l[ e), and 
Algorithm A. 2 'e). We apply the CG method to solve the linear system AI = b. 

1 
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3 
4 
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ro *- b - AI , pq <- r , k <- 0. 
while ||rfe||2 > £ and k < it Max do 

T 

ak < rY^- 

P k Ap k 



h+i <- Ik + ot k p k - 
r k +i 4- Tk - a k Ap k . 

i-fci-fc 

Pk+l <~ Tk+l + PkPk- 

k i- k + 1. 
end while 

I<-h- 
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